#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>

static inline double doublemax(double a, double b) {
    return a>b?a:b;
}

int main(int argc, char *argv[]) {
	double S0, logS, strikePrice;
	double r, realSigma, priceSigma;
	double dt, hedgeInterval, endTime, time, lastHedge;
	double delta, c, d1, d2, timeToMaturity, oldDelta;
	double *dollahZ, *S;
	double mean, variance, difference;
	gsl_rng *randomNumGen;
	int t, timeSteps, run, nruns, var;
	FILE *fp, *fp_commands;
	int i, seed;
	seed = 1337;
	if(argc > 1) seed = atoi(argv[1]);
	fp_commands = fopen("tmp.txt", "w");
	fp = fopen("tmp-data.txt", "w");
	fprintf(fp_commands, "set title \"Hedge interval vs portfolio value\\ntrader sigma = 0.20, stock sigma = 0.20\"\n\
				set xrange [0:32]\n\
				set key right bottom\n\
				set xlabel 'Hedge interval'\n\
				set ylabel 'Portfolio value'\n\
				set output 'hedge_interval.png'\n\
				set terminal png\n\
				plot 'tmp-data.txt' using 1:3 title 'Standard deviation', 'tmp-data.txt' using 1:2 title 'Mean value' \n");
	S = malloc(sizeof(double) * timeSteps);
	dollahZ = malloc(sizeof(double) * timeSteps);
	randomNumGen = gsl_rng_alloc (gsl_rng_mt19937);
	gsl_rng_set(randomNumGen, seed);

	// Input parameters:
	S0 = 100.0;
	strikePrice = 99.0;
	realSigma = 0.20;
	priceSigma = 0.20;
	r = 0.06;
	dt = 0.5; // days
	hedgeInterval = 1.0;
	endTime = 1.0; // years	
	nruns = 500;
	
	
	dt /= 365.;
	
	for(var=1; var<32; var++) {
		hedgeInterval = var*1.0;
		hedgeInterval /= 365.;
		mean = 0.0;
		variance = 0.0;
		for(run = 1; run <= nruns; run++) {
			logS = log(S0);
			timeSteps = round(endTime / dt);
			delta = 0.0;
			dollahZ[0] = 0.0;
			
			timeToMaturity = endTime;
			d1 = (logS - log(strikePrice) + (r + realSigma*realSigma * 0.5) * timeToMaturity)
							/ (realSigma * sqrt(timeToMaturity));
			d2 = d1 - realSigma * sqrt(timeToMaturity);
			c = S0 * gsl_cdf_ugaussian_P(d1) - strikePrice * exp(-r * timeToMaturity) * gsl_cdf_ugaussian_P(d2);

			//printf("Black-Scholes: c = %f\n", c);
			lastHedge = -1.0;
			dollahZ[0] = c;
			for(t = 0; t < timeSteps; t++) {
				time = t*dt;
				S[t] = exp(logS);
				
				// Perform hedging
				//if(1) {
				if(time - lastHedge >= hedgeInterval ) {
					lastHedge = time;
					oldDelta = delta;
					timeToMaturity = endTime - time;
					d1 = (logS - log(strikePrice) + (r + priceSigma*priceSigma * 0.5) * timeToMaturity)
							/ (priceSigma * sqrt(timeToMaturity));
					delta = exp(-r * timeToMaturity) * gsl_cdf_ugaussian_P(d1);
					dollahZ[t] += (oldDelta - delta) * S[t];
				}
				
				// Estimate option value
				d1 = (logS - log(strikePrice) + (r + realSigma*realSigma * 0.5) * timeToMaturity)
						/ (realSigma * sqrt(timeToMaturity));
				d2 = d1 - realSigma * sqrt(timeToMaturity);
				c = S[t] * gsl_cdf_ugaussian_P(d1) - strikePrice * exp(-r * timeToMaturity) * gsl_cdf_ugaussian_P(d2);
				
				//printf("%0.3f, S = %f, $ = %f, delta = %f, port = %f\n", time, S[t], dollahZ[t], delta, dollahZ[t] + delta*S[t] - c);
				//fprintf(fp, "%0.3f\t%f\t%f\n", time, S[t], dollahZ[t] + delta*S[t] - c);
				
				if(t < timeSteps-1) {
					dollahZ[t+1] = dollahZ[t] * exp(r * dt);
					logS += (r - realSigma * realSigma * 0.5) * dt + realSigma * gsl_ran_ugaussian(randomNumGen) * sqrt(dt);
				}
			}
			
			dollahZ[timeSteps-1] += delta*S[timeSteps-1];
			//printf("End delta = %f, selling gives $ = %f\n", delta, dollahZ[timeSteps-1]);
			dollahZ[timeSteps-1] -= doublemax(S[timeSteps-1] - strikePrice, 0.0);
			//printf("After using option, $ = %f\n", dollahZ[timeSteps-1]);
			
			difference = dollahZ[timeSteps-1] - mean;
			mean += difference / run;
			variance += difference * (dollahZ[timeSteps-1] - mean);
		}
		fprintf(fp, "%f %f %f\n", hedgeInterval*365, mean, sqrt(variance/nruns));
	}
	gsl_rng_free(randomNumGen);
	
	return 0;
}
